Towards complete detection of unstable periodic orbits in chaotic systems

We present a rigorous analysis and numerical evidence indicating that a recently developed methodology for detecting unstable periodic orbits is capable of yielding all orbits up to periods limited only by the computer precision. In particular, we argue that an efficient convergence to every periodic orbit can be achieved and the basin of attraction can be made finite and accessible for typical or particularly chosen initial conditions.  2001 Elsevier Science B.V. All rights reserved. PACS: 05.45.Ac; 05.45.Pq Unstable periodic orbits (UPOs) have been recognized as perhaps the most fundamental building blocks of invariant sets in chaotic dynamical systems [1,2]. Many measurable quantities of physical interest can be related to the dynamical properties of the set of infinite number of UPOs embedded in the chaotic set. It is thus of paramount interest in the study of chaotic systems that a complete set of UPOs can be computed. There are several algorithms for computing UPOs in chaotic systems, notably the Biham–Wenzel (BW) method [3] that is applicable to Hénon-like maps [4], the recent Schmelcher–Diakonos (SD) method [5] and its variant developed by two of us (DL) [6]. All these methods appear capable of yielding all UPOs up to reasonably * Corresponding author. E-mail address: [email protected] (Y.-C. Lai). high periods in a variety of model chaotic systems, but so far there is no rigorous assurance that this will be true. The aim of this Letter is to make a step forward by presenting a rigorous analysis of the basin size of the DL method in one dimension and providing numerical evidence in two dimensions suggesting that the method is in fact capable of yielding all UPOs [6]. The key to complete detection of UPOs lies in the basin of attraction for every orbit [7]. That is, there exists a region in the phase space in which initial conditions converge to the UPO under a particular numerical iterative scheme. In order to have a successful detection, this basin must be finite and be accessible for typical (randomly chosen) or particularly chosen initial conditions. Because the DL method is a variant of the SD method, our basin analysis implies that the SD method can also yield a complete set of UPOs. However, as we 0375-9601/01/$ – see front matter  2001 Elsevier Science B.V. All rights reserved. PII: S0375-9601(01) 00 46 37 100 R.L. Davidchack et al. / Physics Letters A 287 (2001) 99–104 will argue heuristically in this Letter, the DL method has the unique property of high efficiency and scalability, which makes a complete detection of UPOs feasible. In a chaotic system, the number Np of UPOs increases exponentially with the period p: Np ∼ eT p , where hT > 0 is the topological entropy of the chaotic set. The traditional Newton–Raphson (NR) method possesses a fast convergent property, but the basin for each individual UPO shrinks exponentially with the period because the number of UPOs increases exponentially with the period. The NR method is thus numerically infeasible for UPOs of high periods [7]. The essence of the SD–DL methodologies lies in a significant enlargement of the basins for UPOs of high periods because they make use of “switches” which make each UPO, in turn, globally convergent. In other words, the SD–DL methods allows for a finite family of iterators, one of which is expected to be globally convergent for each UPO of a given period. On the other hand, the NR method is forced to cope with all UPOs of a given period with the same iterator for all, which thus forces exponentially decreasing basins for each UPO. The convergence of the SD method, however, becomes extremely slow when the period is high [5]. The DL method [6], on the other hand, takes the advantages of, and overcomes the drawbacks of, both the SD and the NR methods so that it is fast and globally convergent. We will argue that, in the DL detection scheme, the basin of attraction to every UPO can be made finite and accessible for random or particularly chosen initial conditions. Numerical evidence also suggests that, starting from all UPOs of the lowest periods, whose detection is guaranteed even by the NR method, the scheme to select initial conditions in the DL algorithm appears to warrant that every UPO of higher periods can be detected. 1 In what follows, we first describe the NR, SD, and DL methods and prove a theorem regarding the finiteness of basin of attraction for one-dimensional maps in the DL scheme. We then argue, with the aid of numerical tracking of basins of convergence, that a similar property can be expected in high dimensions. Since the argument for high dimensions is not rigorous, we provide numerical examples including a hyperbolic chaotic set for which no UPO 1 We still have no rigorous understanding why this is so. is missing, to demonstrate that the DL method can indeed detect all UPOs in an efficient manner. Consider an N -dimensional chaotic map xn+1 = f(xn). The orbit points of period p can be detected as the zeros of the following function: g(x)= f(p)(x)− x, where f(x) is the p times iterated map of f(x). To find zeros of g(x), one usually chooses an initial point x0 and then computes successive corrections xk+1 = xk + δx, which converge to the desired solution. In the NR method [7], the corrections are calculated from a set of N linear equations −J(x)δx = g(x), where J(x)= ∂g/∂x is the Jacobian matrix. The NR method has excellent convergence properties, approximately doubling the number of significant digits upon every iteration, provided that the initial point is within the linear neighborhood of the solution, which shrinks exponentially as the period increases. In the SD method, the corrections are determined as follows: δx = λCg(x), where λ is a small positive number and C is an N ×N matrix, which we call a switching matrix, with elements Cij ∈ {0,±1} such that each row or column contains only one nonzero element. With an appropriate choice of C and a sufficiently small value of λ the above procedure can find any periodic point of a chaotic system. The main advantage of the SD method is that the basin of attraction of each UPO extends far beyond its linear neighborhood, so most initial points converge to a UPO. The iterative scheme in the DL method is as follows: (1) xk+1 = xk + [ 1βg(x)− CJ(x)−1 · Cg(x), where g(x) ≡ ‖g(x)‖ 0 is the length of the vector, and β > 0 is an adjustable parameter. In the vicinity of an UPO, the function g(x) tends to zero and the NR method is restored. Away from the solution and for sufficiently large values of β , the DL scheme is similar to that of the SD method and thus almost completely preserves its global convergence property. Qualitatively, this similarity can be understood by noting that the SD method is in fact the Euler method with step size λ for solving the following system of ODEs: dx/ds = Cg(x), while the DL method is the semi-implicit Euler method [8] with step size h = 1/βg(x) for solving the same system of ODEs. Consequently, with sufficiently small step size, both methods closely follow the solutions to the ODEs and thus share the global convergence property. R.L. Davidchack et al. / Physics Letters A 287 (2001) 99–104 101 Computationally, the DL method for finding all periodic points of period p thus consists of the following steps: (1) Find all periodic orbits of low periods, say, all fixed points and all period-2 orbits, by using any of the described above iteration schemes with a sufficiently large number of randomly chosen seeds; (2) List all 2N ! matrices CN×N , and determine the subset of these necessary for stabilizing arbitrary hyperbolic equilibria; (3) Iterate Eq. (1) by using every periodic point of period p − 1 as seeds, and choosing a matrix C and a number β = β1 > 0; (4) If the sequence {xk} converges to a root of g, then iterate f to find all components of the orbit; (5) Repeat steps (3) and (4) for every matrix C in the list; (6) Repeat steps (3)–(5) for increasing values of β until no more new orbit points are found. A key point of this Letter concerns the scaling of β , in step (6), which leads us to believe with a high degree of numerical confidence that we have a complete set of UPOs, up to a given period. Consider that for very small period-p, once we have saturated β , while we have strong belief that we have collected all the up to period-p points, we can formulate independent verification of completeness using a very large number of randomly chosen initial seeds. Once we have verified our β-saturation threshold for small p, then for higher period-p, by the scalability of our method, that is by following exactly the same steps for each period, we assert with strengthened confidence that our collected list of up-to period-p points is complete. It is demonstrated in Ref. [6] that for two-dimensional maps such as the Ikeda–Hammel–Jones–Moloney (IHJM) map [9], for a high period, say p = 20, the number of UPOs that can be detected as a function of the parameter β saturates at β ≈ 103, indicating that no orbit is missing for β > 103. For higher periods, the value of β needs to be larger. In principle, one can detect, efficiently, all orbits up to periods that are limited only by the computer round-off. Efficiency and the ability to detect complete sets of UPOs are therefore the two features of the DL algorithm. In the sequel, we will provide rigorous analysis and/or numerical evidence for these features. 1. Efficiency. We stress that, with respect to efficiency, the major improved ingredient of the DL method lies in its property of fast convergence, as compared with the SD method. In particular, note that in the neighborhood of a root of g(x), we have ‖g‖ ≈ 0 and, hence, the algorithm is essentially NR whose convergence rate is quadratic, a statement that can be proven rigorously [10]. In addition, in the SD method, vast majority of the computation is spent on finding UPOs that are already found or on trajectories that do not converge at all. In the DL method, the seeds for period-p orbits are the orbit points of all the period(p − 1) orbits, which helps reduce, significantly, the computation time. In actual implementation, we recommend the following seeding scheme: (1) Find all UPOs of low period (say p 3); (2) Choose β ; (3) Use the (p − 1)-period seeds to detect UPOs of period p with all matrices CN×N ; (4) Use these detected period-p UPOs as seeds to detect period-(p+1)UPOs as in step (3); (5) Use the just detected p+ 1 UPOs as seeds for detecting the lower period-p points again; (6) Repeat steps (2)–(5) for increasing β until saturation, i.e., no more new period-p orbits are detected. 2. Complete detection of UPOs. The key requirement for a complete detection of UPOs is that the basin for each UPO not be exponentially small for large p. For a chaotic system, a fairly recent work [7] has shown rigorously that the traditional NR algorithm yields basins whose sizes decrease exponentially with the period. An appealing feature of the DL algorithm is that the basin sizes for distinct periodic orbits can be controlled. In particular, by increasing the value of the parameter β in Eq. (1), the basin can be enlarged. 2 The rest of this Letter is to prove this key statement for one-dimensional systems and provide numerical support for higher-dimensional systems. Theorem. For a one-dimensional map f (x), let g(x) = f p(x)− x and x̄ be a root of g(x). Assume that x̄ is not a critical point, i.e., g′(x̄) = 0. For β 0, define the following function: (2) H(x)≡ x + Cg(x) β|g(x)| −Cg′(x) and also define Nβ,C(x̄) to be the open interval containing x̄ and satisfying |H ′(x)| < 1 for all x ∈ 2 A large value of β, however, requires more time steps for the numerical trajectory to enter the neighborhood of the corresponding periodic orbit. Thus, there is a trade off between enlarging the basins